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Abstract 

The NEMO general circulation ocean model is extended to incorporate three physical 
processes related to ocean surface waves, namely the surface stress (modified by growth 
and dissipation of the oceanic wave field), the turbulent kinetic energy flux from breaking 
waves, and the Stokes-Coriolis force. Experiments are done with NEMO in ocean-only 
(forced) mode and coupled to the ECMWF atmospheric and wave models. Ocean-only 
integrations are forced with fields from the ERA-Interim reanalysis. All three effects 
are noticeable in the extra-tropics, but the sea-state dependent turbulent kinetic energy 
flux yields by far the largest difference. This is partly because the control run has too 
vigorous deep mixing due to an empirical mixing term in NEMO. We investigate the 
relation between this ad hoc mixing and Langmuir turbulence and find that it is much 
more effective than the Langmuir parameterization used in NEMO. The biases in sea 
surface temperature as well as subsurface temperature are reduced, and the total ocean 
heat content exhibits a trend closer to that observed in a recent ocean reanalysis (ORAS4) 
when wave effects are included. Seasonal integrations of the coupled atmosphere-wave- 
ocean model consisting of NEMO, the wave model ECWAM and the atmospheric model 
of ECMWE similarly show that the sea surface temperature biases are greatly reduced 
when the mixing is controlled by the sea state and properly weighted by the thickness 
of the uppermost level of the ocean model. These wave-related physical processes were 
recently implemented in the operational coupled ensemble forecast system of ECMWF. 


1 Introduction 


Surface waves affect the ocean surface boundary layer (OSBL) through a number of processes, 
but perh aps most visibly through bre aking waves which can be seen as whitecaps on the ocean 
surface ( Monahan . 1971 : Wu . 1979[l . These breaking waves enhance the t urbulence in the 
upper part of the ocean significantly ^Craia and Bannen Il994l: CraiaV Il996h . Waves absorb 
energy and momentu r n from the wind fie l d when they grow and in tu r n release it wh e n the y 
break f Janssen et al . 2004 : Rascle et ai . 20061: Ardhuin and Jenkins . 2006 : Janssen . 2012). 
This lowers or raises the stress on the water side (i.e., the stress below the oceanic wave field) 
relative to the air-side stress, depending on whether the sea state is growing or decaying. Only 
when the wave held is in equilibrium with the energy injected by the wind will the stress on 
the two sides of the surface be equal. 

Through the interaction with the Coriolis effect, the Stokes drift velocity associated with the 
wave held adds an additional term to the momentum equation. The ehect was hrst presented 
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bv iHasseJmanm (Il97nh and has since been inves tigat ed for idealized cases bv \ Weher\ (jl983|), 
Jenkins ( 1987h . McWilliams and Restrepo (1999) and McWillia,rn.s a,nd Sullivam ( 2000 ) among 


others. The force is varionsly known as the Stokes-Coriolis force or the Hasselmann force 
depending on whether it is considered to be pnrely an effect of the average Coriolis force acting 
on a particle with a Lagrangian ve locity as given by t he mean cnrren t s and the waves or as 
a tilting of the planetary vorticity nPolton et al\. 120051: xBrostrom, et al\. I201411 . The force does 
not directly modify the total mass transport bnt it will alter the distribntion of momentnm 
over the depth of the Ekman layer ^McWilliams and Restrepd Il999l: \Polton\. 120091) . 


studied in a number of sing 

e-column mixed-layer model experiments ( Craia and Banner . 

1994 


McM 

Jlliams and Restreno. 

1999; 

McWilliam,s and Sullivan. 2000 

Burchard. 

2001 

Ka,ntha a,nd Clavson 

2004; 

Mellon a.nd Blnm.bera 

, 2004 

[ R.a,scle et aJ. 2006: Ardhuin a,nd Jenkins. 

2006 

Hua,na et al. 


2011; 

Ja,nssen\. 20121). Several studies have employed large eddy simulations (LES) to inves 


tigate the impact of Langmnir tnrbnlence in the npper ocean nSkvllinastad and Denboi 11995 


McWilliams et al.\.\1997\:\Teixeira and Belchen\2002:\Polton and Belchen\200Ti\ Grant and Belcher 


2009) . and in some cases even direct nnmerical simulations (DNS) have been employed LSullivan et al. 
2004J ). Most of these studies find that waves do indeed seem to have a rather profound im¬ 
pact on the upper part of the ocean, but there is still considerable disagreement about which 
processes are more important. So far there have been few studies of the wave impact on three- 
dimensional ocean circulation models or fully coupled models of the ocean, the atmosphere and 
the oceanic wave field although the potentia l impact of waves on the c li mate system is re cognized 


nBabamin et all 120091: \ Ca,va,leri et all l2012l: \Fan a,nd GriMesl I2014J) . \Fan et a,l\ (120091) demon¬ 


strated the importance of correctly mod elling momentum and en ergy fluxes from the wave held 
to the ocean under hurricane conditions. Fan and Griffie^{2014) found that the introducti on of 
Langmuir tu r bulen ce following the parameterizations by McWilliams and Sullivan f 2000l) and 


Smvth et al\ fl2002h a s well as the parameterization of mixing by non-breaking waves suggested 


by \Qiao et al\ (120041) signihcantly changed the upper-ocean temperature in long-term coupled 
climate integrations. This latter mixing process appears similar to the mi xing due t o the high 
Reynolds numbers o f the orbital motion of non-breaking waves explored bv lRa^anml (120061) and 


Babanin and Hau^ (120091 ). Using a climate model of intermediate complexity, Wabanin et al\ 


(120091) explored three wave-related mixing processes, namely injection of turbulent kinetic en¬ 
ergy from bre aking waves. Lang i nuir c irculation and the aforementioned mixing by non-breaking 
waves. Like Fa,n a,nd GriMes ( 2014t) they found that all three processe s contributed to the 
mixed layer depth and the temperature of the mixed layer. Similarly, \Hua,na et a,l\ (120111) 
coupled WAVEWATCH III ( Tolman et al . 2002) to a version of the Princeton Ocean Model 
( Blumbera and MelloA 1987 ) and demonstr ated an imp r oved s ummertime temperature prohle 
using the non-breaking parameterization oi lQiao et an (120041) . They found very little direct 
impact by the breaking waves on the temperature. 

These wave-driven processes influence the vertical structure of the temperature and current 
fields in the mixed layer in general, and in the upper few meters in particular. This has 
implications for cou pled models as these processes will affect the feedback between the ocean 
and the atmosphere Llamssen et fli.l . l20131) . However, on shorter time scales and at higher spatial 
resolution it is also clear that these processes will influence the drift of objects and pollutants 
on the sea surface or part i ally o r wholly submerged. This has practical i mpor t ance for oil spill 
modelling (Hackett et al\ 2006 ). search and rescue { Breivik and Allen . 20081: Davidson et al\. 


20091 : iBreivik et all 120131) and dispersion of biological material nRohrs et alV 120141) . 


The NEMO ocean model \Madec and the NEMO teana. 120121) has been coupled to the at¬ 
mospheric model with wave forcing from the wave model as part of the ensemble suite of the 
Integrated Forecast System (IFS) of the European Centre for Medium-Range Weather Fore¬ 
casts (ECMWF) since November 2013 (IFS Cycle 40R1). Here we describe the implementation 
of the three wave effects mentioned above in forced (ocean -only) integrations of NEMO using 
forcing from the ERA-Interim reanalysis { Dee et 201 ll) as well as their implementation in 
a fully coupled atmosphere-wave-ocean seasonal forecast system. 
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The paper is organized as follows. Sec |2] describes the processes that have been implemented 
and lays ont their actual implementation in NEMO. Sec |3] describes the results of long ocean- 
only integrations and co mpares with control runs, observations and the ORAS4 ocean reanalysis 
i Balmaseda et ai . 2013h . Sec 0] describes the coupled atmosphere-wave-ocean coupling used for 
seasonal integrations and compares the results of a control run where no direct coupling exists 
between the wave model and the ocean model to a run where NEM O is force d with stresses, 


turbulent fluxes and Stokes drift from the wave model ECWAM f ECMWF . 2ni3h . Sec [5] 


discusses the results and the dehciencies in the existing model setup. Sec [6] concludes and 
makes suggestions for further work on the investigation of wave effects in ocean-only as well as 
coupled atmosphere-wave-ocean models. 


2 Wave effects in the Ocean Surface Boundary Layer 

Introducing wave forcing in an Eulerian ocean model entails communicating the relevant forcing 
helds from a wave model. We start with a brief presentation of spectral wave models and how 
the two-dimensional wave spectrum relates to the fluxes and helds that have a bearing on the 
ocean surface boundary layer. 


2.1 Fluxes and fields estimated from a spectral wave model 


Third generatio n spec t ral wave model s nHasselma,nn et a,l\ 


1994 : \E.is et all Il999l: \Jansseri 120041: \ Tolman et al 


2002 


Talma,n 


. 1991 

Kom,en et a,l. 

2997: 

Cavaleri et al. 


2007ll solve t he ac tion balance equation (see Eq (1.185) bv I Aomen et al\ (119941) and Eq (2.71) 


bv 1 Janssenl (j200J)) for the wave action density N (a function of the Cartesian co-ordinate x, 
frequency / and direction 6) as follows. 


f d du d 
elk (9x 


A.Tliv 

9x dk) 


SL + S'^, + 5 ',. 


( 1 ) 


Here, the right-hand source terms refer to wind input (in), nonlinear transfer (nl), and dissipa¬ 
tion due to wave breaking (ds), respectively. The dissipation term may include shallow-water 
effects and bottom friction. The gradient in frequency represents shoaling and refraction, and 
CO = O' -|- k ■ u is the absolute frequency as seen by an observer standing still, whereas a is the 
intrinsic frequency as seen by an observer moving with the current. We also note that 


du 

'K 


= 


U, 


( 2 ) 


where Cg is the wave group veloctiy vector. The wave action density is related to the wave 
variance density through N = F/a. In deep water with no current refraction Eq ([T]) reduces to 
the energy balance equation. 


dF 


+ V ■ (cgA) — Sij^ + -|- S'ds, 


(3) 


written here in flux form. Note that the source terms in Eq ([3]) are related to those in the 
action balance equation ([T]) as S' = aS'. 


2.2 The air-side stress modified by surface waves 

The presence of an undulating surface affects the roughness felt by the airflow. The atmospheric 
momentum flux to the oceanic wave held is denoted Tin. It is convenient to dehne an air-side 
friction velocity in relation to the total air-side stress, Ta, as 

uI = tJp^. (4) 
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Here, pa is the surface air density. ICharnocM fjl955h was the hrst to relate the roughness of the 
sea surface to the friction velocity, 


u: 


^0 — «CH- 


9 


(5) 


where acu is known as the Charnock constant. iJanssenI (1l989l. Il99l[) assumed that acn is not 
constant but varies with the sea state. 


acH = 


acH 


yi - T|„/r,’ 


( 6 ) 


where den = 0.006 ^BidloA l2012h and the wave-induced stress. Tin, is related to the wind input 
to the wave field as 

poo 

Pw9 I I S[jidcud0. (7) 

Jo Jo ^ 

Here pw is the water density. A diagnostic s pectral tail proportion al to oj~^ has been applied 
above a cutoff frequency ^Komen et all 1 1994 lECMWFl [20131). The wave-modified drag 
coefficient is then 




K 


log^(10/2;o)’ 

where k = 0.4 is von Karman’s constant. Note that the drag coefficient as dehned here is 
related to the 10-m neutral wind speed, 17 ion- This drag coefficient is computed by ECWAM. 
The Charnock parameter ([6]) is the m ain coupling me chanism between the atmosphere and the 
wave field in ITS, in place since 1998 { Janssen . 2004h . 


2.3 The water-side stress modified by surface waves 

As the wind increases, the wave held responds by hrst growing and storing more momentum. 
In this phase there is a net inhux of momentum to the wave held. Then, as the waves mature 
and the breaking intensihes, the momentum hux from the wave held to the ocean starts to close 
on the hux from the atmosphere to the waves. This is the equilibrium state where dissipation 
matches wind in put, also referr e d to a s fully developed windsea (since the waves c a nnot become 


higher ), see e.g. \Kom.en et an fll994]) : \ World Meteorological Oraanizatiom (jl998|); \Holth,uiisen 


(120071) . Finally, as the wind dies down there will be a net outhux of momentum from the wave 
held, almost all of which will go to the ocean. 

If wind input and dissipation in the wave held were in equilibrium, the air-side stress would 
be equal to the total water-side stress. By water-side stress is meant the stress as seen by the 
Eulerian ocean, i.e., th e momentum hux from the waves. However, most of the time waves 
are not in equilibrium flJan.s.seni. l2012l: \ Janssen et ali. l2013h , giving diherences in air-side and 
water-side stress of the order of 5-10%, with occasional departures much larger in cases where 
the wind suddenly slackens. Likewise, in cases with sudden onset of strong winds the input 
from the wind held will be much larger than the dissipation to the ocean, lowering the water¬ 
side stress to values well below 70% of its normal ratio to the air-side stress. The water-side 
stress thus equals the total atmospheric stress minus the momentum hux absorbed by the wave 
held (positive) minus the momentum injected from breaking waves to the ocean (negative). 
Toe = Ta — Tin “ Tds- Here, the dissipat ion source term is assumed to include all relevant 
dissipative processes. This can be written f ECMW^ 2013h 


p27T poo 

Toc=Ta-pwp/ / —{Sin + Sds)dujd6. 


'0 JO 


CO 


(9) 


The stress from waves is archived as a normalized quantity and is applied as a factor to the 
air-side stress in our implementation in NEMO. 
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2.4 Mixing parameterizations 

The TKE equation with Reynolds averages can be written 


De 

m 


9 


zf / I 11 

= —u'sP' - U[u) 


dui 


P 


^ dxi dxi ^ 


9 /—N 1 
(<e) - 


d 


pwdx, 


Kp') - e- 


( 10 ) 


Here, e = q'^/2 = is the TKE per unit mass (with q the turbulent velocity) and e the 


dissipation rate (see e.g. \Stula fllOSSh . p 152). NEMO has the option of m odelling the evol ution 


of TK E with local closure (a prognostic equation in e only, see l^ftt/ZlIlOSSl pp 203-208 and iPope 


2000l pp 369-373). Assuming that the advective terms are small in comparison and making the 
gradient transport approximation where turbulent coefficients are proportional to the gradients 
in the mean quantities, we arrive at 


— - K 
at 


3 3f 


(11) 


This is the standard one-equation formu lation for NEMO (see the reference manual for NEMO 
v3.4, Madec and the NEMO tea^ 2012 . pp 176-177). Here is the mixing length. The buoy¬ 
ancy term is assumed proportional to the local Brunt-Vaisala frequency. 


N^ = - 


gdp^ 


p dz 

and the shear production is related to the shear of the mean flow. 


( 12 ) 


S" = 


/ (9u 


V dz 


(13) 


Finally, the r nixing len gth is given by a relation hy \Blanke and Deleclusa fll993h . see also 
Caspar et al. (1990) and Madec and the NEMO team f 2012h . pp 177-179. 

Two non-standard mixing processes present in NEMO’s TKE scheme warrant our attention. 
The hrst is an artihcial boost to the TKE known as the ET AU parameterization which is pegge d 
to the surface TKE with an exponential vertical decay (see ldfadec and the NEMO team ( 2012 1. 
Sec 10.1), 

er{z) = 0.05ei exp 2 ;/hT-. (14) 

The depth scale h-r can vary with longitude from 0.5 m at the Equator to 30 m poleward of 44° 
or be hxed at 10 m. The coefficient, here 0.05, can also be varied. The sec ond rn i xing p rocess 
of interest to us is a parameterization of Langmuir turbulence according to AxelH (2002) which 
has been implemented in NEMO. The vertical velocity wlc of tho Langmuir cells is assumed 
to peak at i7Lc/2, half the maximum depth to which Langmuir cells penetrate. 


wlc { z ) = CLc^sSin 


nz 


H, 


LC 


, 0 > z > 77] 


LC- 


(15) 


Here Vs is the surface Stokes drift speed and clc = 0.15 is a coeffic ient. \AxeJA (1200211. by 


making an analogy with the characteristic convective velocity scale W’Alessio et all 1199811 


further assumed the Langmuir production term in the TKE equation flTU]) could be written 


Pi^ciz) = 


w 


LC 


77] 


LC 


(16) 


Th is production term will at tain a maximum value in the interior of the mixed layer. 


Craia and Bannen (1199411 . CB94 hereafter, demonstrated that as waves break they will con¬ 


siderably mod ify the verti cal dissipation prohle from the traditional law-of-the-wal l where dissi¬ 


pation (X z ^ LStuli 1198811 . With wave breaking CB94 found dissipation cx z \ Terra,v et al. 
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fll996h and \Drennan et an (119961) later demonstrated that the observed dissipation rates nn- 
der breaking waves are indeed much higher than anticipated by the la w of the wall. C B94’s 
model has since been extended to a two-equation turbulence model by Burchar^ ( 2001 ). who 
demonstrated that the injection of turbulent kinetic energy from breaking waves was sufficient 
to successfully model the evolution of the mixed layer representative of North Sea conditions. 
CB94 suggested that the flux of turbulence kinetic energy (TKE) should be related to the water 
friction velocity tc* as 

*hoc = PwacB^^^^ (17) 


CB94 assumed that ocb was a constant ~ 100, but n oted that its range would p robably be 
between 50 and 150, depending on the sea stat e (see also Mellor and Blum,hero 20041) . The TKE 
flux f rom breaking waves is related (see e.g. Janssen et al. 2004 Rascle et al. 20061: Janssen 


2OI2I : I Janssen et a/.ll2013f) to the dissipation source function of a spectral wave model as 


r-27r 


*hoc = -pw9 


Sj, dud0 = -p^mul 


(18) 


For consistency we have written the energy flux from the waves (thus always negative), m ~ 

— ■\/pa./ PwttcB, normalized by the air friction velocity n*. In NEMO the ene rgy flux from break¬ 
i ng wa ves is introduced as a Dirichlet boundary condition on TKE, following lMeJor a,nd Blumhera 
(120041) ■ It is assumed that in the wave-affected layer the mixing length can be set to a con¬ 
stant /w = where the surface roughness length relates to the significant wave height as 
= 0.5iJs, and that in this near-surface region diffusion balances dissipation. In this case th e 
TKE equation takes a simple exponential solution [see Eq (10) bv lMe//or a,nd Blumhera\\200^ . 


e(z) = eoexp(2Az/3). 


Here the inverse length scale is 


A = l3/(SgBK^)]^^^z. ^ 


(19) 


( 20 ) 


with Sg = 0.2 and B = 16.6 given by Mellor a,nd Yamada ( 1982h . This is how the flux from 
breaking waves is implemented in NEMO v3.4. This allows the following simple boundary 
condition 

^ -.2/3 l^ocl 


Co — - (15.8 q;cb)" 


( 21 ) 


However, the inverse depth scale (I2UD is sea state dependent, and for a wave height of, say, 
2.5 m, which is close to the global mean, ~ 0.5 m. Thus, e(z) varies rapidly with depth, and 
we have modified the boundary condition (l?ID by weighting the surface value by the thickness 
of the topmost level to attain an average value more representative of the turbulence near the 
surface of the model, 

1 

Cl = — / e(z) dz. (22) 


'-L 


Here L = Azi/2 is the depth of the T-point of the first level. Integrating Eq (IT^ is straight¬ 
forward, and the average TKE boundary condition (l22|l becomes 


Cl — Co 


2XL 


[1 — exp(—2AL/3)]. 


(23) 


It is clear that as the vertical resolution increases the difference between ei and Cq becomes 
smaller, and in the limit the two coincide. The weighting (I23p is thus less important with 
higher vertical resolu t ion. It is worth noting that the exponential profile (IT^ assumed by 
Mellor and Blvmhera (120041) is only valid very near the surface, and in fact CB94 had already 
found the solution to the more general case where the mixing length is allowed to vary with 
depth. We have not implemented this operationally, but preliminary tests suggest that the 
effect is to roughly double the depth over which the TKE from breaking waves is distributed. 
The derivation is presented in the appendix. 
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2.5 The Stokes-Coriolis forcing 


Waves set up a Lagrangian displacement Vg in the down-wave direction known as the Stokes 
drift velocity ^StokesL 118471 ) . Although it decays rapidly with depth, it can be substantial near 
the surface (|vs| ~ 0.7ms“^). In combination with the Earth’s rotation it adds an additional 


veering to the upper-ocean currents known as the Stokes-Coriolis force ^Hasselmannl 119701) . 

1 dr 


Du 1 N 

-— = —Vp + u + Vg X /z + - —, 
Dt p p oz 


(24) 


Here / is the Coriolis frequency, z is the upward unit vector, p is the pressure and t is the 
stress. The full two-dimensional spe c trum is in prin ciple required to compute the Stokes drift 
velocity profile f Janssen et al . 2004 : Janssenl. 2012), 


-27r 


VgU 


= Air 


/ke2‘'F(/,»)d/d«. 


(25) 


This is computationall y demanding an d full two-dimensional wave spectra from a numerical 


wave model (see e.g. \ECMWF] l2013l) may not always available. It i s therefore cus t omar y 


to introduce a simp l ihed, monoch r omat ic Stokes drift p r ohle (see e.g., \Carniel et al\ (120051) : 
Polton et fli.l (120051) : \Saetra et al\ (1200711 : xTamura et al\ (1201211 1. However, it was shown by 
Breivik et al\ (120141 ) that this profile is a poor match to the full profile and that the following 
parameterization gives a considerable improvement. 


-. 2 keZ 


\e[Z) = Vo 


1 — SkpZ 


(26) 


Here the subscript “e” distinguishes the approximate prohle from the full Stokes drift velocity 
prohle (12^ . The surfa ce Stokes drift ve locity vector vq is computed by ECWAM and is available 
both in ERA-Interim ( Bee et ai . 2011 1 and from the operational ECMWF forecasts ( ECMWF . 

201 oi l. 

To compute the prohle (l26ll we must hnd the i nverse depth scale kp. This is r elated to the 
transport Tg throu gh the exponential i ntegral Ei nAbramowitz and Steauru . Il972l) and can be 
solved analytically ( Breivik et al . 20141) to yield 


T = 

-‘-s — 


|vo|eV"Ei(l/4) 

8kp 


Rearranging we get the following expression for the inverse depth scale, 

|vo|eV4Ei(l/4) 


kp = 


Here Ei{l/A) ^ 1.34, thus 


kp 


8Tg 


|vo| 

5.97Tg 


The n-th order spectral moment is dehned as 


^27r /‘CO 




/"^(/.Ojd/dO. 


(27) 


(28) 


(29) 


(30) 


The mean freque ncy is dehned as / = mi/mo n World Meteorological OraanizatiorL Il998l: 
Holthuiisem . 120071 ) and the signihcant wave height Big = Ay/mo- We can derive the hrst moment 
from the integrated parameters of a wave model or from wave observations and hnd an estimate 
for the Stokes transport, 

„ 27r^ 


16 




(31) 


Here kg = (sin 6'g, cos 6'g) is the unit vector in the direction Og of the Stokes transport. W e 
approximate the Stokes transport direction by the surface Stokes drift (see \Breivik et a/.l (1201411 ). 
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3 Ocean-only Forced Model Experiments 


The NEMO model is run on a tripolar ORCA 1° grid configuration with 42 vertical levels. The 
uppermost level is 10 m thick. The model is c o upled to LIM2, a two - Rvel thermodynamic- 


dynamic sea ice model ( Fichefet and Maqueda, . 1997t Bouillon et al. . 2009) and is relaxed 


weakly towards a climatology in temperature (3-yr e-folding time). No sea surface temperature 
(SST) relaxation is performed. The ORCA grid is such that the re solution is increased towards 
the E quator (roughly 1/3°) to better resolve tropical waves (see \Madec and the NEMO team 
( 2012 ) for details on the ORCA grid). 

The atmospheric and wave forcing fie l ds ha ve been computed from the ERA-Interim re¬ 


analysis \Simmons et all 120071: \Dee et all 1201111 . ERA-Interim is a continuously updated at 


mospheric and wave field reanalysis starting in 1979. The resolution of the wave model is 1.0° 
on the Equator but the resolution is kept approximately constant globally through the use 
of a quasi- regular ktitud e-longitude grid where grid points are progressively removed toward 
the poles Uansseru. 1200411 . Similarly, the atmospheric model fields are archived on a reduced 
Gaussian grid of approximately 0.75° resolution at the Equator. Some care has to be taken 
when interpolating between these grids, in particular where wave parameters are interpolated 
from the ECWAM grid to the ORCA grid. NEMO requires fluxes to be defined in all ocean 
points. However, there are discrepancies between the ice coverage and the land-sea m ask of the 
wave grid and the ocean grid. This is solved by reverting to the ECMWF drag law Uanssem. 
2008; Edson et al . 2013h where ECWAM has ice or land. 


C'D(;^ = 10m) = (a + 6f/f/)/f/fo^ 


(32) 


The coefficients are a = 1.03 x 10“^, b = 0.04 x 10“^, pi = 1.48 and p 2 = 0.21. Here Uio is the 
10-m wind speed from ERA-Interim. Where ECWAM and NEMO agree on open water, the 
stress is computed from the drag coefficient of ECWAM, 


Ta — PaCjOwU] 


ION- 


(33) 


Here, t/ioN is the neutral 10-m wind speed, available on the ECWAM grid. The conversion to 
water-side stress is implemented as 

Toe = TTa, (34) 

where f is the ratio of water-side to air-side stress (see also Eq IMjl . It is this parameter which 
is archived by ERA-Interim. 

A standard integration period covering the ERA-Interim period from 1979 up until the end of 
2009 has been used in the following. A summary of the settings for the model runs can be found 
in Tables [T] and [21 Four experiments with the new wave-related effects are presented, all com¬ 
pared against a control experiment, CTRL, where standard settings are used for NEMO. The 
CTRL experiment includes a parameterization of the TK E flux from breaking waves ( CB94 ) 
without explicit sea state information (see Eq flT7[l and also lMadec and the NEMO team fl2012fl . 
Sec 10.1) but has no averaging over the topmost model level fl22D . The stress in CTRL is com¬ 
puted using the ECMWF drag law (1521) with air-side stress (1551) . In all runs the NEMO VAR 
observation operator y = 'H(x) relating model state x to observation space (y) is applied. 
This allows a comparison of model integrations (in our case without assimilation) against the 
large number of q uality-controlled tempera t ure a nd temperature-salinity profiles compiled in 
the EN3 data set ( Inalebv and HuddlestoA . 2007 1. Seasonal averages (December to January, 
DJF, and June to August, JJA) over the period 1989-2008 are used to compare SST fields in 
the following. 

The first wave experiment, TAUOC, uses the water-side stress described in Sec [23] together 
with the ECWAM drag coefficient (see Sec 12.2[) . The effect is confined mostly to areas with 
rapidly developing weather systems in the extra-tropics (FigHj), where the sea state will be quite 
far from equilibrium. There is a slight weakening of the wind stress along the west coast of 
South America and along the coast of south-west Africa, leading to decreased upwelling. This 







































is mainly a consequence of differences between the ECWAM drag coefficient and the drag law 
fl3^ caused by limited fetch near the coast. The overall effect is a slight reduction of the bias 
in the tropics compared to EN3 near-surface temperature measurements (0.1 K, not shown). 
In the tropics, differences are most likely due to differences in the drag over swell compared to 
the drag law (1521) used for the CTRL experiment. 

The second experiment, WTKE, introduces the TKE flux dTS]) from ECWAM (described 
in Sec 12.4j) . The differences found in the extra-tropics amount to more than 2 K (Fig [2]). The 
large difference suggests that the standard settings of NEMO overdoes the mixing due to waves, 
especially with coarse vertical resolution (cf Eq [22]) . 

The STCOR experiment introduces Stokes-Coriolis forcing from ECWAM as described in 
Sec 12.51 The largest impact (Fig[3]) is found in areas with extra-tropical cyclones in combination 
with strong temperature gradients, such as across the Gulf Stream and the Kuro-Shio. 

Another experiment called LOW where a Law-of-the-wall boundary condition is applied, 
i.e. no TKE flux from breaking waves, was also performed. The motivation is to establish a 
lower bound on the mixing. The difference between the WTKE and LOW runs is much smaller 
(not shown) than the difference between CTRL and WTKE. 

We now combine the three wave experiments TAUOC, WTKE and STCOR in one exper¬ 
iment referred to as the WAVE run. Fig H] shows the standard deviation (upper curves) and 
the biases of the CTRL run (blue) (blue) (blue) (blue) (blue) (blue) (blue) (blue) and a wave 
run including all three wave effects (green). The most striking feature is the large-amplitude 
annual cycle observed for the bias of the CTRL run. This amplitude is much weaker for the 
run with wave effects, and the seasonality is not at all so clear. Due to the huge differences to 
the mixing in the runs with and without wave effects, the heat uptake of the ocean also differs 
signihcantly. Fig [5] shows the global total heat content anomaly relative to 1979 (the start of 
the period) for the experiment with wave effects (green) v the CTRL experiment (blue). We 
also show the ocean reanalysis OR AS4 (red) which is based on an earlier version of NEMO 
(v3.0), see \B almaseda et al\ (120131). with observat ions assimilated using the NEMOVAR 3D- 
VAR assimilation system (iMooen.sen et an. I2ni2al). The r eanal ysis covers the period 1957 to 
present, with forcing provided by ERA-40 ( Uppala et all 2005 ) and later ERA-Interim from 
1979 and onwards. The similarity with ORAS4 is clear, and the trends almost identical. It 
is also clear that the impact of the wave mixing establishes itself within the first two years of 
the integration. The pronounced annual cycle (dashed lines) is due to the difference in oceanic 
volume in the southern and northern hemispheres. 

To further a ssess the impact on th e surface temperature we now compare with OLv2, an 
SST analysis by Reynolds et al. ( 20021) . The comparison with the CTRL run in Fig [5^ reveals 
large biases, especially in the summer hemisphere. These biases are reduced when wave effects 
are included (Fig[6)D), especially in the northern extra-tropics. This is mainly due to a more 
correct level of mixing. The long-term SST fields from ORAS4 are very similar to OI_v2 since a 
strong relaxation to the O I_v2 gridded SST was ap plied as a flux correction between December 
1981 and December 2009 nBalmaseda et al\. 120131 ) and are not shown here. 

Finally, to test the relative impact of the ETAU TKE boost (I14p and the Langmuir tur¬ 
bulence parameterization (1T6|) . we switched off the Langmuir circulation (NOLC) and ETAU 
(NOETAU) in two additional experiments. As is evident from Fig[71 the Langmuir turbulence 
has only modest impact on the temperature in the mixed layer in the extra-tropics (50 m depth 
in the northern extra-tropics is shown in Fig [7]). Switching off ETAU has a much larger impact 
on the mixing, both in terms of bias and random error. 


4 The Coupled Atmosphere-Wave-Ocean Model 

Coupling between the wave model and the ocean model was hrst implemented operationally 
in the ensemble suite of the IFS in Cycle 40R1 in November 2013, but was limited to the 
mixing (WTKE) and Stokes-Coriolis forcing (STCOR). The seasonal integrations described 
here include in addition the modihed stress (TAUOC) and the ice model LIM2 but are run at 


9 

































lower spatial resolution than the operational ensemble forecasts. Fields are exchanged every 
three hours (the coupling timestep) between IFS/ECWAM and NEMO. The atmospheric model 
is run with a spectral trnncation of T255 (corresponding to ronghly 78 km) with 91 vertical 
levels to 1 Pa. ECWAM is rnn at 1.5° resolution while NEMO is run on the ORCA 1° grid 
described in SecO Fig [H] shows a flow chart outlining the sequence of execntion of the various 
components of the IFS-ECWAM-NEMO single executable over two coupling timesteps, which 
can be snmmarized as follows. 


1. The atmospheric component (IFS) is integrated (internal time step 2,700 s) one conpling 
time step (10,800 s), yielding wind fields for ECWAM as well as radiation and evaporation 
minus precipitation fields for NEMO 

2. ECWAM (internal time step 900 s) is tightly conpled to IFS and returns sea snrface 
ronghness to the atmospheric boundary layer at every atmospheric time step. After one 
NEMO conpling time step, ECWAM also gives Stokes drift velocity, tnrbulent kinetic 
energy and waterside stress to NEMO 

3. NEMO is integrated (internal time step 3,600 s) one conpling time step, yielding SST and 
snrface cnrrents to IFS 

4. All model components have now been integrated one conpling time step and the sequence 
begins anew 


Surface currents and SST are commnnicated back to the atmospheric model and will affect the 
stress and the temperature of the boundary layer. This in tnrn affects the oceanic wave field, 
bnt there is presently no direct feedback from NE MO to ECWAM. The co npling between the 
different components is described in more detail by Moaensen et al. f 2012bl) . Here we compare 
a setnp for seasonal integrations to seven months, with three ensemble members starting from 1 
May and 1 November for the period 1993-2012. The CTRL experiment is rnn with a standard 
CB94 type wave mixing which corresponds to acs = 100 (Eq|2T]), similar to the ocean-only 
CTRL experiment presented in Sec [3l The wind stress is compnted using the ECMWF drag 
law fl5^ and 10-m wind vectors from IFS. The wave experiment includes the three processes 
TAUOC, WTKE and STCOR described in Sec [2] and Sec El Fig [9] reveals large differences in 
the bias in the northern extra-tropics relative to ERA-Interim for the boreal snmmer (JJA) at 
a lead time of one to three months from 1 May. Panel (a) shows the bias of the CTRL rnn, 
with large cold biases in the northern extra-tropics. These biases are broadly similar to what 
was fonnd from the ocean-only (forced) runs presented in Sec El (cf FigE]). Panel (b) reveals 
the biases in the run with wave effects to be generally much smaller, although there is a certain 
deterioration of the npwelling area along the coast of Baja California. Note that the cold bias 
in the eastern eqnatorial Pacihc (the “cold tongue”) is reduced slightly. The results are similar 
for the southern hemisphere summer (DJF), although the bias reduction is not as strong as for 
the northern hemisphere. The seasonal variation of the bias fonnd for the forced (ocean only) 
runs in Fig Elis also present in the coupled integrations, see Fig [101 Again, the bias is greatly 
reduced in the wave run. 


5 Discussion 

We have introduced three wave effects in NEMO, namely the sea-state dependent water-side 
stress, the energy flux from breaking waves and the Stokes-Coriolis force. Using ocean-only 
integrations and experiments with a coupled system consisting of the atmospheric model IFS, 
the wave model ECWAM and NEMO, we demonstrated that the impact of the wave effects 
is particularly noticeable in the extra-tropics. Of the three processes, the modification of the 
mixing (WTKE) has the largest impact (FigE]), but as we discuss below, this is also related 
to the additional mixing found in NEMO. The impact of the modified stress (TAUOC) and 
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Stokes-Coriolis (STCOR) is also significant (on the order of 0.5 K locally, see Figs [U and 
|3]). It is also important to note that compared to the law-of-the-wall experiment (LOW), the 
WTKE differs by only 0.5 K. This again snggests that the CTRL experiment has too vigorous 
mixing. In ocean-only integrations we see a reduction of the temperature bias in the mixed 
layer, particularly in the extra-tropical summer (FigHj). This manifests itself in a more realistic 
oceanic heat uptake (Fig|5]). The coupled seasonal integrations show a similar reduction in 
bias (compared to ERA-Interim, Fig [9]) as the ocean-only wave run (see Fig fTOjl . The mixing 
is strongly influenced by the ETAU parameterization fflTll . It is clear that the TKE scheme 
fllip has too shallow mixing without this parameterization, and it is also clear that the present 
parameterization of Langmuir turbulence flT^ does very little due to its vertical structure which 
has the distinguishing feature that it will put most if not all of the enhanced turbulence deep 
into the mixed layer and nothing near the surface. It is thus unable to transport down enough 
heat to make a substantial difference. This explains why its impact on the temperature in the 
OSBL (Fig [7]) is so much smaller than that from the ETAU term. Its vertical prohle ffT5]l is 
very dif ferent from Langmuir pa r ameterizations involving the shear of the Stokes drift velocity 
profile ( McWilliams et al. . 1997 : Polton and Belcher . 2007 : Grant and Belcher . 2009), 




<9Vs 

dz 


(35) 


Here, is the horizontal velocity fluctuations and w' the vertical. Due to the strong shear 
of the Stokes drift profile this would add a larger contribution near the surface. The ETAU 
profile (1TT|1 is quite similar to fl35|l and it appears to act as a parameterization for Langmuir 
turbulence with its characterist i c exponential decay with depth (l3^, o r similarly mixing by 
non-breaking waves »Qiao et all l2004 iBabanim 120061: iHuana et a/.l. 1201111 . It facilitates deeper 
penetr ation of mixing from surface processes than w hat is normally assumed from breaking 
waves A Orant a,nd Belchen. 120091: (Belcher et all l2012h . 


6 Concluding remarks and further work 

The ocean-only integrations and coupled seasonal integrations all suggest that the right level 
of mixing is very important for reducing the temperature bias in the upper part of the ocean 
and also for the oceanic heat uptake. An important result is that introducing wave-enhanced 
mixing must be done in such a way that the thickness of the uppermost layer is accounted for. 
This is done with the present implementation by weighting with the thickness of the layer and is 
essential with model configurations with a thick uppermost level, e.g. ORCA1L42 as discussed 
here. This would not be necessary if a flux boundary condition had been used for the TKE from 
breaking waves, but this is not the case in NE MO, which uses the surface boundary condition 
hrst proposed bv \Mellor a,nd Blumbercn (120041) . The impact of mixing on SST is clearly shown 
in Fig [6] where runs with and without wave effects are compared with the OI_v2 SST analysis. 
That the seasonal cycle of the CTRL run is distorted by too vigorous mixing is clear, and FigH] 
shows that the annual cycle in biases extends well below the surface. The conclusions from the 
ocean-only experiments that temperature biases are reduced by introduction of wave-induced 
mixing are borne out by the seasonal coupled integrations which essentially show the same 
bias reduction (Fig [9]). It is important to note here, though, that the additional ad hoc deep 
mixing in NEMO interacts with the surface processes and that without this additional mixing 
the model fails to mix deeply enough. We speculate that this mechanism is really masking 
Langmuir turbulence or mixing from non-breaking waves. More work is clearly needed with 
ocean circulation models and coupled models to fully answer the question of which mixing 
processes are dominant in the OSBL, but it is clear that getting the mixing right is a balancing 
act between the right deep mixing and the right mixing near the surface, and these processes 
are probably all wave-related. 

These results are relevant for assessing the impact surface waves have on climate projections 


nBabanin et all 120091: iFan and GriffieA\2Q14i] . and a natural next step would be to investigate 
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the impact of waves on long, decadal to centnry-wide inte grations fsee a lso th e Co-ordinated 
Ocean-Wave Climate Projections (COWCLIP) initiative, Hem^er et al. f 2ni2[p . One candi¬ 
date for forcing ocean-on ly integrations would be the recently completed ERA 20th cen tury 
reanalysis (ERA -20C, see Poli et aZI 2013 . Hershach et al. 2013 . de Boisseson et al. 2014 , and 
Dee et oZI 2014). This opens up the possibility of running ce ntury-long N EMO integrations 
with wave effects from a state-of-the-art version of ECWAM ^Bidlou 120121) since all relevant 
parameters have been archived in the new reanalysis. For coupl ed climate projections, the 
nearest candidate would be EC-Earth ( Hazeleaer et al . 2010l 2012) which operates a modihed 
version of an earlier cycle of IFS. Such experiments would help determining the importance of 
waves in the climate system, rather than just the impact of climate change on the wave climate. 


A The dissipation profile 

The exponential prohle flT^ for the balance of 


^ fl Q 

dzV'^^’^dz) Bl 


(36) 


assumed by lMe//or and Blumberal (120041) is only valid very near the surface where the mixing 
length can be assumed constant = kz^. CB94 presented the solution to the more general case 
where the mixing length is allowed to vary with depth. This equation has a power-law solution 
[cf CB94, Eq (23)], 

e{z) = eo (37) 


z^ — z 


where 


n = 


-j 


1/2 


S„tBB J 


= 2.4. 


This leads to a slightly more complicated expression for the vertical average ([22]), 

-2n/3+l' 


Cl — Co 


L(2n/3 - 1) 


i + T 


(38) 


(39) 


than Eq (l23|) . The consequence is that the mixing pe netrates about twice as dee p as in the 
case where the exponential approximation assumed bv IMellor and Blvmbem (120041) . 
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Horizontal grid 

ORCA 1° 

Vertical resolution 

42 levs (10 m top lev) 

Time step 

3600 s 

Time period 

1979-2009 

Atmospheric forcing 

ERA-Interim 

Data assimilation 

OFF 

SST damping 

OFF 

3D damping to dim 

ON (3 yr Newtonian relaxation) 

Bulk parameterization 

COARE 

Ice model 

LIM2 


Table 1: Overview of the settings common to all forced (ocean only) experiments. 
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Physical process/parameterization 

Experiment description 

Stress 

TKE flux 

Stokes-Coriolis 

Langmuir 

ETAU 

CTRL: 

Control experiment 

Drag law 

Eq ([22D 

CB94 TKE flux 
ttCB = 100 Eq fITTD 

Off 

On 

On 

TAUOC: Water-side 

stress 

ECWAM 
stress (ED 

as CTRL 

Off 

On 

On 

WTKE: Sea-state 
dependent TKE flux 

as CTRL 

ECWAM 

TKE flux fllSp 

Off 

On 

On 

STCOR: Stokes-Coriolis 
forcing 

as CTRL 

as CTRL 

ECWAM 
Stokes (I2TD 

On 

On 

WAVE: 

as TAUOC 

as WTKE 

as STCOR 

On 

On 

All three wave effects 






LOW: Law-of-the-wall 

as CTRL 

Off 

Off 

On 

On 

NOLC: Langmuir off 

as CTRL 

as CTRL 

Off 

Off 

On 

NOETAU: ETAU off 

as CTRL 

as CTRL 

Off 

On 

Off 


Table 2: Overview of the settings of the ocean-only (forced) experiments. Departures from the 
CTRL experiment marked with bold. 
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Figure 1: Long-term SST differences between a run with water-side stress modulated by the 
ECWAM wave model (TAUOC experiment, see Table 2) and the CTRL run. 
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Figure 2: Long-term SST differences between a run with TKE flux from ECWAM (WTKE 
experiment, see Table 2) and the CTRL run. Note that the color scale is ±2 K. 
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Figure 3: Long-term SST differences between a run with Stokes-Coriolis forcing (STCOR ex¬ 
periment, see Table 2) and the CTRL run. 
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Figure 4: A comparison of su rface (a) and 50 m depth (b) EN3 temperature observations 
ninalebv and Huddlestom. 120071) in the northern extra-tropics in the CTRL run (blue) and a 
run with all three wave effects switched on (green), see Table 2 for more information about the 
experiments. The upper curves show the standard deviation while the lower curves represent 
the bias. A 90-day running mean is employed. 
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Figure 5: Panel a: The global ocean heat content anomaly (relative to the start of the time 
series) in the upper 300 m. The run with all three wave effects switched on is shown in green 
(see Table 2). Full lines represent a 12-month moving average. The ORAS4 reanalysis is shown 
in red and the CTRL run in blue. Panel b: Total ocean heat content anomaly. 
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Figure 6: Panel a: Long-term SST differences between the W.evnolds et al\ (1200211 OI_v2 SST 
analysis and the CTRL run. Panel b: Differences between the run with wave effects and the 
CTRL run. 
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Figure 7: A comparison of the impact of ETAU, the parame terizat i on for enhanced deep mixing 
in NEMO, and the Langmuir mixing parameterization bv \Axela (120021) . The northern extra¬ 
tropics (dehned as north of 20° N) at 50 m depth is shown here. The blue curve shows the CTRL 
run where both processes are switched on (default). Switching off the Langmuir mixing (red) 
is seen to have a much smaller impact than switching off the ETAU parameterization (green). 
Finally, a run where both processes are switched off (black) shows that the combined effect is 
again dominated by the ETAU parameterization. Mode led temperature is compared to EN3 
temperature observations ( Inalehv and Huddleston . 2007|l . The upper curves show the standard 
deviation while the lower curves represent the bias. A 90-day running mean is employed. 
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Time 

Figure 8: A flow chart of the coupling between the components of the single executable IFS- 
ECWAM-NEMO model setup. Two coupling time steps are shown to illustrate the sequence. 
This is how the operational ensemble prediction system has been run since Cycle 40R1 (Novem¬ 
ber 2013). The seasonal integration experiments described here contain in addition the LIM2 
ice model and are run at lower atmospheric resolution (T255). First the atmospheric component 
(IFS) is integrated (internal time step 2,700 s) one coupling time step (10,800 s), yielding wind 
helds for ECWAM as well as radiation and evaporation minus precipitation helds for NEMO. 
ECWAM (internal time step 900 s) is tightly coupled to IFS and gives sea surface roughness 
for the atmospheric boundary layer at every atmospheric time step. After one NEMO coupling 
time step, ECWAM also gives Stokes drift velocity, turbulent kinetic energy and waterside stress 
to NEMO. NEMO is integrated (internal time step 3,600 s) one coupling time step, yielding 
SST and surface currents to IFS. All model components have now been integrated one coupling 
time step (10,800 s), and the sequence begins anew with the next coupling time step. 
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Figure 9: Panel a: SST bias relative to ERA-Interim for the coupled seasonal CTRL run (JJA). 
Panel b: Same for the wave run. 
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Figur e 10: A comparison of surface EN3 temperature observations Unalebv and Huddleston 


20071 ) in the northern extra-tropics for coupled integrations of IFS-ECWAM-NEMO as a func¬ 
tion of lead time. The start date is 1 May. A three-member ensemble is run to a lead time of 
seven months for the years 1993-2008. The dashed red (lower) curve shows the bias of the run 
with wave effects, and the green dashed curve shows the bias of the CTRL run. The upper full 
lines show the root-mean-square (RMS) error. 
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